

NASA Technical Memorandum 78759 



- 4 ' 



* 




f 


SPINEQ: A Program for Determining 
Aircraft Equilibrium Spin 
Characteristics Including Stability 



William M. Adams, Jr. 


NOVEMBER 1978 







NASA Technical Memorandum 78759 


SPINEQ: A Program for Determining 
Aircraft Equilibrium Spin 
Characteristics Including Stability 


William M. Adams, Jr. 
Langley Research Center 
Hampton, Virginia 


IWNSA 

National Aeronautics 
and Space Administration 

Scientific and Technicai 
information Office 


1978 


CONTENTS 


SUMMARY 1 

INTRODUCTION 1 

SYMBOLS 1 

MATHEMATICAL MODEL AND SOLUTION TECHNIQUE 4 

Coordinate Systems and Transformations 4 

Mathematical Characterization of an Equxlibrium Spin 6 

System Linearized About an Equilibrium Spin Condition 10 

PROGRAM DESCRIPTION 11 

Overview of Program Flow 11 

Functional Description of Major Subroutines and 

Definition of Key Parameters 12 

SPINEQ Input Requirements 21 

SPINEQ Output Identification 23 

PROGRAM CHANGES FOR NEW CONFIGURATIONS 28 

SPINEQ (Main Program) 29 

AERO 29 

STABLE 29 

POTENTIAL SPINEQ APPLICATIONS AND MODIFICATIONS 29 

PROGRAM AVAILABILITY 30 

APPENDIX - INPUT AND OUTPUT FOR A SAMPLE CASE 31 

REFERENCES 82 


lii 


1 


SUMMARY 


A computer program, SPINEQ, is described which can algebraically solve the 
nonlinear equations of motion for equilibrium spin conditions. Linear charac- 
teristics of the airplane about the equilibrium points are also determined. 

The theoretical basis of the program is outlined, computational flow is 
shown, the functions of major subroutines are described, and key parameters 
directing the computations are identified. Program input and output are 
described and illustrated by means of a test case. The program is available 
from COSMIC. 


INTRODUCTION 

This paper describes a computer program that was written to determine equi- 
librium spin characteristics. It is a modified version of the program that was 
employed to obtain the results presented in reference 1. Prior to the method 
described in reference 1, the analytical methods described in the literature 
for determining airplane equilibrium spin characteristics were the numerical 
integration of nonlinear equations of motion until approximately steady condi- 
tions were achieved (refs. 2, 3, 4, and 5) and the solution of a simplified set 
of equations for approximate equilibrium spin conditions (refs. 6, 7, and 8) . 

In reference 1, a technique for determining equilibrium spin conditions by uti- 
lizing nonlinear programing techniques was described which, for a rigid airplane 
and given nonlinear aerodynamic data, could determine precise equilibrium spin 
conditions. Stability characteristics were also determined by linearization of 
the equations of motion about an equilibrium spin condition. 

Included herein are brief descriptions of the mathematical model and the 
solution technique; a flow chart of the program; a description of the function 
of major subroutines, including key parameters controlling program operation; 
a description of input and output; and a test case, including sample input and 
output . 


SYMBOLS 

A,B matrices in linearized state equation ^ = A^ + Bu 

b wing span 

c mean aerodynamic chord 

F external force vector 

aerodynamic force vector 


f 

G 


g 


g 


h 

I 


J 


J 


n 


J 


nn 


J 


K 






Mo 


iti 


p,q,r 


/\ /s A 

R,T,Zj 


R 

S 


Ti^T2 


vector of nonlinear functions yielding f = (d,3/ V/V,p,q,r) 
vector of nonlinear functions yielding = (d,3,V/VQ,p,q,r,0,^) 
acceleration due to gravity (assumed constant) 
constraint function 

upper and lower limits, respectively, on constrained variables 
gradient of constraint functions with respect to H 
maximum (g - g^, g^ - g) 
aircraft altitude 

inertia tensor with moments and products of inertia computed about 
center of mass 

function to be minimized in order to find equilibrium spin condition 

vector of partial derivatives of J with respect to vector q 

matrix of second partial derivatives with respect to q 

- T 

augmented function, J = J + g Kg 

diagonal matrix of weights that multiply constraint violations 
(described in discussion of DAVIDN subroutine) 

rolling moment about X body axis 

pitching moment about Y body axis 

yawing moment about Z body axis 

mass of airplane 

angular rates about body axes 

1 

dynamic pressure, 

A 

unit vectors in cylindrical coordinates; Zj is directed toward 
center of Earth (fig. 3) 

radius of helical path of airplane 

wing area 

matrices displaying coordinate system transformations (see 
figs. 1 and 2) 
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t time 

U control vector, u'^ = 

u deviation from nominal control vector, U(t) " Uf^{t) 

u,v,w components of velocity vector (see fig- 2) 

V airspeed 

X,Y,Z axes 

x,y,z position triple (when devoid of subscripts, a body-axis coordinate 
system is referred to) 

a angle of attack 

B angle of sideslip 

Y flight-path angle, tan 

^e'^^a'^r elevator, aileron, and rudder deflections (positive 6g is trailing 
edge down, positive 6^ is right trailing edge down, positive 6^. 
IS trailing edge left) 

= (d,<i),r,R,’p,ip') 

= (a,B,V/Vo,p,q,r,9,(t)) 

5 _ deviation from nominal value of 5 

p atmospheric density 

\p,6,(p angles defining transformation between inertial and body axes 

(see fig. 1) 

\p' angle from horizontal projection of airplane longitudinal axis to 

horizontal component of V, positive counterclockwise when looking 
up (see fig. 3) 

I 2 2 2 

w angular rate about center of mass, \Jp + q + r 

Subscripts : 

H horizontal component 

I referred to inertial coordinate system 

N nominal value 

o evaluated at initial time 
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s stall angle 

w referred to wind axis system 

Special symbols: 

{"') unit vector 

("^) vector 

{■^) tensor of second rank 

(•) derivative with respect to time 

( ) transpose of matrix 

( ) * value which minimizes J 

{ } column vector 

MATHEMATICAL MODEL AND SOLUTION TECHNIQUE 

The airplane equations of motion and an algebraic method of solving them 
for equilibrium spin conditions by using nonlinear programing techniques are 
discussed in this section. The solutions that are obtained satisfy the require- 
ments of a steady developed spin if the equilibrium points are stable. Solu- 
tions having linear representations that are unstable either correspond to 
oscillatory spin conditions or are not actual developed spin conditions. 


Coordinate Systems and Transformations 

Three coordinate systems are employed in this analysis. These are an iner- 
tial set, a set fixed in the body, and a wind axis set. The rotations involved 
in transforming from one set of axes to another are shown in figures 1 and 2. 



Yaw Pitch Roll 

Figure 1.- Inertial axes to body axes transformation (numbered 
subscripts refer to intermediate axes) . 
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u 



Figure 2.- Wind axes to body axes transformation. 

In figure 1, the airplane is viewed toward the origin from a point on a positive 
axis. Transformations are as follows: 


(a) Inertial axes to body axes (z^- is in the direction of the force due 
to gravity^ 






Tl = 


cos Ip cos 0 

cos \p sin 0 sin (p 
- sin ^ cos cj) 

cos Ip sin 0 cos c}) 
+ sin \p sin (J) 


sin ip cos 9 

cos ip cos (j) 

+ sin ip sin 0 sin (p 

sin Ip sin 0 cos (p 
- cos Ip sin (j) 


-sin 0 
cos 0 sin (j) 

cos 0 cos (j) 


(b) Body axes to wind axes 



X 


= T2(a,3) \ y / 


where a (angle of attack) and 3 (angle of sideslip) are wind incidence 
angles positive as shown in figure 2 and 


^2 = 


cos 

3 

sin 

3 

“1 

0 


cos 

a 

0 

sin 

a 

-sin 

3 

cos 

3 

0 


0 


1 

0 


_ 0 


0 


i_ 


-sin 

a 

0 

cos 

a_ 
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Mathematical Characterization of an Equilibrium Spin 

Equilibrium spin conditions can be stated mathematically as follows . 
Figure 3 will aid in the interpretation of some of the terms . 



Figure 3.- Top view of helical path. 


a > ttg 


-y • /\ 
0) = IpZj 


(Angle of attack is above stall) 


Angular velocity of airplane about its\ 
center of mass is vertical / 


( 1 ) 

( 2 ) 


V = Ri|;T - hZj /Tra 3 ectory xs along a helical path with\ 

V vertical helix axis (t = Zj. x r) / 


(3) 



= (») 


(constant angular velocity) 


(4) 


->• /s 

F = -mRU) R 


/Force IS such that helical motion\ 
\ persists / 


(5) 


In these equations, thrust effects have been ignored and variation of density 
with altitude has been neglected. 
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The aerodynamic forces in SPINEQ are expressed relative to a body-fixed 
nonprincipal set of axes for which the equations of motion are given by 

i ->■ 

mV + coxmV = F = F^+ mgZj (6) 

I*a) + o)xli.a)) = M (7) 


Here F is composed of forces of aerodynamic and gravitational origin. The 
moments M are due to aerodynamic forces only. The aerodynamic force coeffi- 
cients and their derivatives are defined in the tabulations that follow: 










Pitching moment 


C 


m 


My 




3 c 


m 


36. 



3 Cm 


3 


2 E. 

2V 


X body axis 


Z body axis 

force component 


force component 

_ 


c - 

L 7 

3Cx 

■ 86 e 


c - 

3<5e 


These coefficients are assumed to be nonlinear functions of a (angle of 
attack) and, in some cases, of B (angle of sideslip) . A precise definition 
of the form assumed for these coefficients is given as follows. Instructions 
are given in a later section for modification of this form when necessary. 

~ Cjj(a,S/U=0) + Cx^ (ot,3)5e 

Cy = Cy(a,B/U=0) + Cy~ (a,3)*5g + Cy. (a,B)<5g 

Oe ua 

+ Cy^ (a, 3) 6^ + — |Cy^(a)p + Cy^(a)^ 


Cz = Cz(a,3,u=0) + Czg^(a,3)6e 

Cl = c^(a,3,u=0) + Cl (a,3)6g + c. (a,3)6g 

+ Cj^^(a,3)6^ + |^[CZp(a)p + Ci^(a)rJ 

Cm = c„,(a,3,u=0) + Ci„^^(a,3)6e + c^^Wq 
c^ = Cj^(a,3,u=0) + Cj,^^(a,3)6g + c^^ (a,3)5a 
+ Cng (a,3)i5r + ^^n^(a)p + c^^(a)^ 
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In SPINEQ, the force equations have been transformed to the wind axis 
system and solved for 6t, 0, and V. A search is then made for the aircraft 

state which meets the helical motion requirements of equations (2) and (3) and 
the following conditions for equilibrium: 



The number of independent variables required to specify this state is reduced 
from 12 to 6 by satisfying the equality constraints of equations (2) and (3) 
explicitly and by arbitrarily setting ip(t*) =0. A search is made over the 
set of variables 



where R and ip' are defined in figure 3 and Y is the flight-path angle. 
The search, therefore, is for a vector q* for which 
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The search for 


ri* xs made by fxnding the h = r)* which mxnimizes 


J(T1) 


D 

= ^ f^2 = (5^2 + g2 ^ 


V2 


+ p2 + q2 + r-2 


i=l 


svibject to the previously mentxoned constraints. If a > and J(r|*) = 0, 
an equilibrium spin has been found. 

The solution for ri* is obtained in an iterative manner, beginning with 
an initial estimate The method of search is initially equivalent to a 

gradient procedure, but second-order information is accumulated during the 
iterations and quadratic convergence is approached near the end of the itera- 
tions. Thus, the method initially has the sureness of convergence toward a 
local minimum of a gradient procedure. Near the solution, the rapidness of 
convergence of a second-order method is approached without computing and invert- 
ing the matrix of the second partial derivatives at each point. The method 
employed is Davidon's (ref. 9) , as modified by Fletcher and Powell (ref. 10) . 

A refinement phase is then entered which uses the accumulated second-order 
information to perform a search by using Newton's method. The program which 
served as the basis for SPINEQ is described in reference 11. 

The user may occasionally encounter convergence difficulties. A local 
minimum may be found for which J(n*) ^ 0, in which case an alternate initial 
guess should be employed. The program may also fail to converge as a result of 
discontinuities in slope which arise from the interpolation routine chosen to 
find the aerodynamic coefficients. Three options for performing the interpola- 
tion are provided. One of these is a cubic spline technique which produces con- 
tinuous first derivatives. 


System Linearized About an Equilibrium Spin Condition 

A linearized representation of the airplane equations of motion is devel- 
oped in order to obtain some information about the stability of the equilibrium 
spin conditions. The linear representation neglects density variation with 
altitude and only includes those states which effect the forces and moments. 
Define 5 such that 


"'^(t) = (a,g,V/VQ,p,q,r,e,^) 


where 


z(t) = G[E(t) ,U(t)] 
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to first order 


kit) . =(t) - =„(t) - (v=g)J 5 + (5u=)^&(t) - U„(t)] 

or 

5{t) = A(t)^(t) + B(t)u(t) 

The elements of A and B are shown explicitly in appendix C of reference 1. 
For an equilibrium condition, A and B are constant matrices. 

All stable eigenvalues indicate a return to the equilibrium spin conditions 
for sufficiently small deviations. Additional analysis would be required to 
determine whether a solution with an unstable linear representation is an 
oscillatory or a no-spin solution. 

PROGRAM DESCRIPTION 
Overview of Program Flow 

Mathematical statement .- The mathematical statement of the problem solved 
is as follows; 

Object: Find the q* which minimizes 

j(n) n is N X 1 

subject to the constraints 

g(n) = {o} g is m X 1 

where m is the number of elements of g. 

Approach : Form 

m 

j = J + ^ K^jmax 
i=l 

Find the q = q* which minimizes J as for i = l,m. The q* which 

minimizes J(q) as solves the problem (ref. 12). 
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Numerical implementation .- The computations performed in obtaining a solu- 
tion are indicated in figure 4. This figure presents an overview of the compu- 
tational flow. Some subroutines that are contained in the program are not 
shown . 


Functional Description 
and Definition of 


of Ma]or Subroutines 
Key Parameters 


SPINEQ (Main Program) 

Purpose : 

(1) Directs computation 

(2) Input of all data here (i.e., aerodynamic tables, mass, inertia, geo- 

metric characteristics, weights, tolerances, and boundaries associ- 
ated with constraints and initial solution estimate) 

(3) Decides upon the number of DAVIDN phases and whether to refine the 

solution by using Newton's method 

(4) Calls STABLE to get linear characteristics 
Subroutines called: 

Information obtained 
I Minimum of augmented function 

V" “ (^nn)'" 

Full output of information of interest 
at end of DAVIDN phase 

Solution for = {o}, g = {o}; treat- 
ing all violated constraints as 
equality constraints 

Density at altitude specified by input 
(ref. 13) 

Linear representation of system about an 
equilibrium point 

I Eigenvalues and eigenvectors of linear. 

V system 


Name 

DAVIDN 

FFXGGX 

NEWT 

AT62 

STABLE 
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Figure 4.- SPINEQ flow chart. 






























Key parameters: 


Name 

N,MEQ,MINEQ 

ISOFA = 1 
15 

IALPHA,IBETA 

DELF 

IPART 


IPERM 

NUMB 

DAVIDN 
Purpose : 


Descriptxon 

Number of xndependent varxables, equalxty 
constraints, and xnequalxty constraxnts, 
respectively 

See discussion of FFXGGX 

Dummy vari^le used to decide when to end 
run (xf 15 = 999 execution is ended) 

Type of interpolatxon to perform in deter- 
minxng the aerodynamic coeffxcxents 

Convergence criterion for DAVIDN 

Evaluatxon of gradient and full output 

Evaluation of gradxent and reduced 
output 

No gradient evaluatxon and reduced 
output 

Reorderxng of g array (violatxons 
moved to top) 

No reorderxng of g (some reorderxng 
could have been done prevxously) 

Number of cases that have been studied 
so far in a given run 


IPART 




r 

1 
0 

r 


r. 


^0 


IPERM = 


l>0 


(1) Finds the minxmiam of the augmented function J, given a set of wexghts 

associated wxth the constraxnt violatxons 

(2) Determxnes an estxmate of the xnverse of the second derivative of the 

augmented function J 



(3) Outputs descrxptxve parameters 
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Subroutines called: 


Name 


INT 


CUBIC 


MATIN 


Information obtained 

Value of augmented function at particular 
points 

Augmented gradient at start of each one- 
^ dimensional minimization 

^Coefficients are determined which fit a 
cubic through last four points evalu- 
ated during a one-dimensional search 

Point (lying between the data points) at 
^ which this function has a minimum 

Determinant of H to check that it 
remains positive definite. This is 
an expensive bit of information which 
could, perhaps, be deleted except when 
convergence difficulties are occurring 


An overview of DAVIDN operation: 

(1) An initial step size for the one-dimensional search is selected (CCAY) . 

(2) A direction of search is defined. 

(3) A one-dimensional search is performed which finds the scalar a such 

that J(n - Oi HJ^) is the minimum of J in the direction of 

(4) The metric H is augmented. (See ref. 10.) 

(5) Several stopping criteria are examined: 

(a) If the number of one-dimensional minimizations performed is less 

than N, the dimension of r|f continue the iterations. 

(b) If the number of one-dimensional minimizations performed is 

greater than 120, stop (failure to converge) . 

(c) If |j(r) - a HJp) - J(r|) | < 10“^^, a minimimi has been 

found. 

(d) If I [j(ri - a HJ^) - J(r|)j / J(n) j < DELF, a minimum has 

been found. 

Other stopping criteria are suggested in reference 10 and could easily be 
added . 
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(6) If another one-dimensional minxmization is required, the initial step 

IS taken to be equal to the a found in the previous one-dimensional 
search and steps (2) to (5) are repeated. 


Key parameters: 


Name 


N 

CCAY 


Description 

Number of variables in n vector 

Initial step size taken at the start of 
each one-dimensional minimization 


ICNT 

ITT 

DELF 

CAYMIN 

P 

I PART 
IPERM 

CUBIC 

Purpose : 


Number of one-dimensional minimizations 
performed 

Total number of function (J) evalua- 
tions made during a particular Davidon 
phase 

Stopping criteria referring to relative 
change in J between succeeding one- 
dimensional minimizations 

Scalar which defines the step, in a 

specified direction, which minimizes J 

Below statement 2003, this is a vector 
in the direction of one-dimensional 
' search 

Below statement 31 in DAVIDN, P is 
^ actual step made 

Defined previously in description of 
SPINEQ 

Defined previously in description of 
SPINEQ (>0 is not allowed because it 
would prevent the effect of a new con- 
straint violation from entering) 


(1) Fits a cubic equation to data generated during the one-dimensional 

search 

(2) Finds the value of a, the scalar search parameter, for which this 

cubic has a minimum in the vicinity of the four data points used to 
construct the cubic 
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Subroutxne called: 


Name Informatxon obtained 

MATIN Coeffxcxents of cubxc equatxon 

MATIN 

Purpose : 

(1) Fxnds the xnverse of a square matrxx and xts determxnant 

(2) Solves a matrxx equation Ax = y 

INT 

Purpose : 

(1) Evaluatxon of the augmented functxon J 

(2) Evaluatxon of the gradxent of the augmented functxon xf a one- 

dxmensxonal search is bexng xnxtxated 

Informatxon obtained 
Value of unaugmented functxon 
Constraxnt vxolations 
Gradxent of unaugmented function 
Gradxent of constraxnt vxolations 

Descriptxon 

Number of elements of r) 

Niamber of constraxnt violations 
Previously defxned xn SPINEQ 

FFXGGX 
Purpose : 

(1) Determxnes the values of the constraxnt and unaugmented function values 


Subroutxne called: 
Name 

FFXGGX 

Key parameters: 
Name 
N 

MC 

I PART 
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(2) Reorders the elements of the constraint vector when required to move 

those that are violated to the top portion of the constraint vector 
array 

(3) Computes, when required, the gradient of the unaugmented function and 

of the constraints 

Subroutine called: 

Name Information obtained 


SETUP 


Key parameters: 


Vector of constraints g 

Vector f used to compute unaugmented 
function value 


Description 


I SOFA 


IPERM 


I SOFA = 


IPERM = 


1 Increment is defined in elements 

of n for perturbing for genera- 
tion of numerical partials 

2 Increment has been defined; not done 

again 

<0 All constraint violations are moved 
to top of g array 

>0 g is not reordered from what it was 
before this entry into FFXGGX. It 
may have been reordered previously. 


Purpose : 


Takes the vector r) and determines f and g. 


Subroutine called: 


Information obtained 


Aerodynamic forces and moments 
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Key parameters: 


Name 


Description 


n 


Vector of variables 


Inertial and geometric 
data 


P 

lu'^1 


AERO 


Purpose : 


See list of symbols 


Atmospheric density 

Upper and lower limits, respectively, on 
constrained variables 


(1) Interpolates aerodynamic tables for aerodynamic force and moment 
coefficients 


(2) Computes aerodynamic forces and moments 
Subroutines called; 


Name 


Information obtained 


lUNl 


CSIUNI 


IBI 


CSIBI 


Interpolated aerodynamic coefficients 
that are a function of a only for 
linear or quadratic fit 

Interpolated aerodynamic coefficients 
that are a function of a only for 
cvibic spline fit 

Interpolated aerodynamic coefficients 
that are a function of a and 3 for 
linear or quadratic fit 

Interpolated aerodynamic coefficients 
that are a function of a and 3 for 
cubic spline fits 
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Key parameters: 


Name Description 

IALPHA,IBETA Variables defining the type of 

interpolation : 

Linear 

Quadratic 

Cubic spline 

NEWT 


IALPHA,IBETA = \2 
3 


Purpose : 

Determine variables and Lagrange multipliers at which the first derivative 
of the augmented function vanishes and at which the constraint violations 
(all taken to be equality constraints) vanish. 

Subroutines called: 


Name 

MATIN 


FFXGGX 


Key parameters: 
Name 
MC 


IPERM 


I PART 


Information obtained 

- 

Inverse of gHg 

Value of unaugmented function 

Values of gradient of unaugmented 
function 

Values of selected constraint variables 
Values of gradient of these violations 


Description 

Number of violated constraints (remains 
constant throughout the Newton-Raphson 
phase (0 < MC < N) 

Must be >0, which means that if a new con- 
straint violation occurs during itera- 
tions, it will not be considered 

Defined in SPINEQ; -1 is not allowed 
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STABLE 


Purpose : 

(1) Gets linearized representation of the system about an equilibrium point 

(2) Gets eigenvalues and eigenvectors of the resulting system matrix 


Subroutines called: 


Information obtained 


Aerodynamic forces and moments 


Aerodynamic force and moment coefficients 


Eigenvalues and eigenvectors of system 
matrix 


SPINEQ Input Requirements 


Variables 


Format 


Description 


IALPHA,IBETA,IAEROUT,IDEG 415 IALPHA,IBETA Type of fit with respect 

to a and B 


1 Linear 


IALPHA,IBETA = (2 Quadratic 

Q3 Cubic spline 


lAEROUT = 


0 No aerodynamic output 

1 Aerodynamic tables output 


IDEG = 


0 ALPHAT,BETAT are input in 

radians 

1 Convert ALPHAT,BETAT to 

radians 


ALPHAT(I) , 
I=1,NPTSA 


9F8.5 


Tabular array of a at which 
aerodynamic data are given 


BETAT(I) , 
I=1,NPTSB 


9F8.5 


Tabular array of B at which 
aerodynamic data are given 


{ ( (CD(I,J,K) ,I=1,NPTSA) , 9F8.5 

J=1,NPTSB) ,K=1,NTABD) 


NTABD tables of aerodynamic 
coefficients that are func- 
tion of both a and B 
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Variables 


Format 


Description 


( (CS (I, J) ,I=1,NPTSA) , 9F8.5 NTABS tables of aerodynamic 

J=1,NTABS) coefficients that are func- 

tions of a 

NOTE: See the listing of subroutine AERO to find the necessary correspondence 

between CD and CS and the physically meaningful coefficients 
(e.g., Cj^jCyfCg) and arrange your tables to correspond. 

B,ELL,S,EM 4E20.8 Span, mean aerodynamic chord, 

wing area or appropriate 
reference area, total mass 

EIX,EIY,EIZ,EIXZ 4E20.8- ^XX'^YY'^ZZ'^XZ the 

center of mass 


GRAV,C0NFAC1, 4E20.8 

C0NFAC2 


(SU(I) ,1=1,9) 4E16.8 


(SL(I) ,1=1,9) 4E16.8 


GRAV Acceleration of gravity 

CONFACl Factor to convert Hq to feet 
for AT62, default is 1 

C0NFAC2 Factor to convert density in 
slugs/ft^ output by AT62 
to desired units, default 
is 1 


Upper limits on variables 
which may be constrained 


SU{1) 

~ ^max 

SU(6) 

= V 
max 

SU(2) 


SU(7) 

“ ^e ,max 

SU(3) 

^max 

SU(8) 

~ '^a,max 

SU(4) 

^max 

SU(9) 

_ r 

,max 

SU{5) 

= ip' 
^max 



Lower 

limits 

on the 

previously 


defined variables which may 
be constrained 


NOTE: In SPINEQ runs to date, constraints have generally been placed only 

upon Y- The user may wish to redefine the constraint variables, in 
which case changes would be required in both SPINEQ and SETUP. This 
would very likely be the case if the program were modified to compute 
trim conditions for maneuvers. 
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Variables 


Format 


Description 


(CKAY(I) ,1=1,9) 4E16.8 Weights corresponding to 

previous constraint 
possibilities 

(GEPSS (I) , 1=1 , 9) 4E16.8 Tolerances on violations of 

constraints 

N,MEQ,MINEQ 315 Number of variables, equality 

constraints, and inequality 
constraints, respectively 

NOTE: No provision has been made in SPINEQ for equality constraints so MEQ 

must be zero. If it is desired in a modification to include equality 
constraints, changes would be required in SPINEQ and SETUP. 

HO E20.8 Altitude at which equilibrium 

spin conditions are desired 

15 15 Stopping parameter at the 

beginning of a loop for mul- 
tiple cases. If 15 ^ 999, 
another case is to be run 

(XBARO(I) ,1=1,6) 3E20.9 Initial estimate for equilib- 

rium spin conditions 

DE,DA,DR 3E20.8 Fixed control settings for 

elevator, ailerons, and 
rudder, respectively 

An example of the input for a test case is presented in the appendix. The 

core required for this run was 105 000 octal to load and 72 000 octal for ]ob 

execution. Job execution time was 18 sec on a Control Data CYBER 173 Computer 

System. Job execution cost was two dollars. 


SPINEQ Output Identification 

The outputs from each subroutine are listed in the order that they occur 
within that subroutine. This does not, in all cases, correspond to the actual 
order of output. 
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SPINEQ 


(1) lALPHAjIBETA Types of aerodynamic fit: 

Linear 
Quadratic 
Cubic spline 

Convert a and 3 to radians 

a and 3 have been input in 
radians 

Output aerodynamic tables 
No aerodynamic table outputs 


If lAEROUT = 0, skip (3) . 

(3) ALPHAT ARRAY 
BETAT ARRAY 
FUNCTIONS of a a 
FUNCTIONS of 3 

(4) B,L,S,MASS 

(5) IX,IY,IZ,IXZ 

(6) GAMMA, ALPHA, OMEGA, 
BETA,PSIPR,V 

(7) DEMAX, DAMAX,DRMAX 

(8) DEMIN, DAMIN,DRMIN 

(9) H0,GRAV 

(10) C0NFAC1,C0NFAC2 


Output a and 3 values where 
aerodynamic data are defined and 
3 then output the aerodynamic data 

Span, mean aerodynamic chord, wing 
area, and mass, respectively 

^XX'^YY'^ZZ'^XZ' respectively, 
about center of mass 

Upper and lower constraints on y, 
a, 0 ), 3, , and V, 

respectively 

Upper limits on 6g, 6^, and 6^,, 

respectively 

Lower limits on 6^, 6^, and 6^,, 

respectively 

Altitude and gravitational 
acceleration, respectively 

Conversion factors; defined in 
list of input; related to input 
to and resulting output from 
AT62 


( 2 ) IDEG , lAEROUT 


IALPHA,IBETA = / 2 


IDEG = 


lAEROUT = 
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(11) 

DENSITY 

(12) 

N,MEQ,MINEQ 

(13) 

THETA , PHI , GAMMA , 
RADIUS, PSIDT,PSI 

(14) 

DE,DA,DR 

(15) 

GEAR, TOLERANCE , WEIGHT 

(16) 

FXX 

FFXGGX 


(1) 

(GBAR(I) ,I=1,MC) 

(2) 

(FX(I) ,I=1,N) 

SETUP 

(No output unless IPART 

(1) 

RA , PS IDT , GAMMAH , THETA 

(2) 

PHI,VH,W,V 

(3) 

P,Q,R 

(4) 

ALPHAD,BETAD 

(5) 

E1,E2,E3 

(6) 

E4,E5,E6 

DAVIDN 


(1) 

CYCLE 

(2) 

(ALPHA (I) ,I=1,N) 


Density at which an equilibrium 
spin condition is sought 

Number of variables, equality con- 
straints, and inequality con- 
straints, respectively 

Initial estimate for equilibrium 
spin conditions (6,(j),Y,R,'i',>p' ) 

Fixed control setting for 
elevator, aileron, and 
rudder, respectively 

Constraints, tolerances, and 
weights, respectively 

H matrix at the end of a DAVIDN 
phase 


Output of constraints (only output 
if IPART > 0) 

Gradient of unaugmented function 
(only output if IPART = 0) 


R,^,Y>0 

(}),Vh,W,V 

P/q^r 

a and 3 (in degrees) 

d,3,V/V 

P,q,r 


Number of one-dimensional itera- 
tions performed in this Davidon 
phase 

Current estimate of the n which 
minimizes the function 
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(3) DFDK 


(4) (FX(I) ,I=1,N) 

(5) KMIN 

(6) FI 

(7) F2 

(8) DELXFX 


(9) DETERM 


(10) IT 


Inner product of J^| and vector 
defining the direction of search 

Gradient of the augmented function 

Scalar at which J (l) + KMIN*P) 
is a minimum in direction of vec- 
tor P (p is not normalized) 

Value of J before current one- 
dimensional minimization 

Value of J at new p after one- 
dimensional minimization has 
been performed 

Inner product of previous direc- 
tion of search and gradient 
evaluated at point where one- 
dimensional minimum was found. 
Product should be zero if deter- 
mination of the one-dimensional 
minimum were precise 

f 

The first time, this is the (1,1) 
element of the H matrix, which 
should be positive 

The second time, it is the deter- 
minant of H which should be 
^ positive 

Number of function evaluations 
called for to date in DAVIDN 


(11) H 


Matrix which approximates the 

inverse of second partial deriv- 
atives of augmented function 


NEWT 

If MC = 0, skip down to output (12). The variable MC is the number of 
equality constraints plus violated inequality constraints at the end of the 
DAVIDN phase. 

(1) AMBDA Initial estimate of the 

Lagrange multiplier. 

AMBDA (I) = KBAR(I) *GBAR(I) 

(2) FALPHA Gradient of J with respect to rj 
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(3) ALPHA 


(4) GEAR 

(5) ICOUNT 

(6) D 

(7) ALPHA 

(8) F 


Initial estimate of r|* (the mini- 
mizing n upon entry into NEWT) 

The vector of violated constraints, 
all of which are treated as 
equality constraints 

Counter of the number of itera- 
tions performed in NEWT 

Product of ( g H g^)(g H g'^) ; 

computed to check on goodness 
of inverse 

Current estimate of n* 


(9) AMDA 


Current estimate of Lagrange 
multipliers 


(10) FALPHA 


J 


n 


(11) GEAR 


Current values for constraint 
violations 


Sequence (7) to (11) is repeated until convergence is achieved. One skips 
from the beginning to output (12) if MC = 0. 


(12) 

FALPHA 

Initial value of 

(13) 

ALPHA 

Initial estimate of q* upon entry 
into NEWT 

(14) 

ICOUNT 

Counter of the number of itera- 
tions within NEWT 

(15) 

ALPHA 

Current estimate of q* 

(16) 

F 

Current value of J 

(17) 

FALPHA 


Sequence (14) to (17) 
STAELE 

is repeated until convergence is achieved. 

(1) 

CONTROLS 

'^e'‘^a''^r 

(2) 

THETA, PH I, GAMMA 

Qf<P,Y 
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(3) RADIUS, PSIDT,PSI PRIME 

(4) ALPHA, BETA, V 

(5) P,Q,R 

(6) x,y,z 

(7) TRANSPOSE OF THE GRADIENT 
OF F WITH RESPECT TO 
CONTROLS 

(8) TRANSPOSE OF THE GRADIENT 
OF F WITH RESPECT TO 
STATE 

(9) ROOTS OF FX 

(10) ROOTS IN POLAR FORM 

(11) EIGENVECTORS IN POLAR 
FORM 

AERO 

(1) ALPHA, BETA 


(2) ALPHAT 


(3) BETAT 


(4) COEFF 


a,6,v 

P/q^r 

Coordinates of R 



Eigenvalues of F^ 

Eigenvalues of F^ expressed in 
polar form 

Eigenvectors expressed in polar 
form 


a and S values at initial call 
to AERO 

a points at which data are avail- 
able. This IS output only for 
the initial call to AERO 

3 points at which data are avail- 
able. This is output only for 
the initial call to AERO 

Values of the aerodynamic coeffi- 
cients at the initial call to 
AERO 


Outputs 1 to 4 are debugging aids and could be deleted whenever desired. 
Output for a test case configuration is presented in the appendix. 


PROGRAM CHANGES FOR NEW CONFIGURATIONS 

When a new configuration is to be analyzed, some changes will typically be 
necessary in the program. The changes arise because of differences in the com- 
plexity of the aerodynamic model for a new configuration. An effort has been 
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made to Ixmit the resulting programing changes to a few specific locations. 
These locations and the types of changes required will now be described. 


SPINEQ (Main Program) 

Specific values, which are dimensions required for the aerodynamic data, 
are assigned to several variables. No changes in dimension or common state- 
ments are required. 

NTABD Number of coefficients that are a 

function of both a and B 

NTABS Number of coefficients that are only a 

function of a 

NPTSA Number of values of a at which aero- 

dynamic data are available 

NPTSB Niamber of values of B at which aero- 

dynamic data are available 


AERO 

Currently it is assumed that there are 25 aerodynamic coefficients that 
define the aerodynamic representation. A specific correspondence is made in 
AERO between these coefficients and the elements of an array given the name 
COEFF. It might well be that only 20 coefficients will be available for a 
new configuration. The user will have to make the proper correspondence between 
variables such as C 2 and 0^ and his COEFF array and then zero out the 

physical coefficients that are either not available or are neglected. 

If new coefficients other than those now considered are to be added, the 
new coefficients will have to be defined in AERO and passed to STABLE through 
COMMON/CFUN/. 


STABLE 

If new coefficients other than those now considered are to be added, the 
new coefficients will have to be defined in AERO and passed to STABLE through 
COMMON/CFUN/ . 


POTENTIAL SPINEQ APPLICATIONS AND MODIFICATIONS 

SPINEQ can provide rapid, accurate solutions of the nonlinear airplane 
equations of motion for equilibrium spin conditions. In addition, the stabil- 
ity characteristics of a linear representation of the aircraft are determined 
about an equilibrium condition. Of course, the realism of the solutions is 
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dependent upon the accuracy of the aerodynamic data. The program can be uti- 
lized to study the sensitivity of equilibrium spin conditions to variations in 
selected aerodynamic or inertial data. It can also be used to aid in evalu- 
ating the adequacy of aerodynamic data in the spin region when accompanying 
experimental spin tests are available. 

The aerodynamic data available vary from configuration to configuration 
insofar as the number of coefficients measured are concerned as well as their 
functional dependence upon angles of attack and sideslip. For example, for one 
configuration, a coefficient may be a nonlinear function of a and 6 whereas 
for a new configuration, it may be assumed to be a function of a only. For 
these reasons, a few minor changes will typically be required each time a new 
configuration is to be analyzed. The required changes are discussed in an 
earlier section entitled "Program Changes for New Configurations." 

SPINEQ could be employed to compute equilibrium conditions for steady 
maneuvers. This would require the addition of thrust and control deflections 
as design variables. Generalization of the relationships between the body axis 
angular rates and the Euler angle angular rates would be required. Alternate 
choices for the design variable set and the set of constrained functions would 
also be necessary. 


PROGRAM AVAILABILITY 
The program is available through COSMIC. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
August 28, 1978 
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INPUT AND OUTPUT FOR A SAMPLE CASE 
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